import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 设置支持中文的字体，这里以'SimHei'为例。如果你使用其他字体，请替换font_name中的字体名称。
font_name = 'SimHei'  # 或者使用你系统中的其他中文字体名

# 更新matplotlib的字体属性，使得其支持中文显示
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = [font_name]

# 如果需要支持负号(-)的正常显示，可以添加以下配置
plt.rcParams['axes.unicode_minus'] = False

# ================= 参数设置 =================
S0 = 100      # 初始股价
K = 100       # 执行价
r = 0.05      # 无风险利率
sigma = 0.2   # 波动率
T = 1         # 到期时间（年）
N_paths = 1000   # 每次模拟的路径数
N_reps = 500     # 重复模拟的次数（用于生成分布）

# Black-Scholes 理论价格（用于参考）
def black_scholes_call(S0, K, r, sigma, T):
    d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    C = S0 * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)
    return C

C_BS = black_scholes_call(S0, K, r, sigma, T)

# 存储每次重复模拟的估计值
C_MC_list = []  # 标准蒙特卡洛估计
C_CV_list = []  # 控制变量法估计

np.random.seed(42)  # 保证结果可复现

for _ in range(N_reps):
    # 生成随机数
    Z = np.random.standard_normal(N_paths)
    
    # 模拟到期股价
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
    
    # 到期收益
    payoff = np.maximum(ST - K, 0)
    
    # 蒙特卡洛估计
    C_MC = np.exp(-r * T) * np.mean(payoff)
    # 在列表里保存本次模拟结果
    C_MC_list.append(C_MC)
    
    # 控制变量法
    E_ST = S0 * np.exp(r * T)  # S_T 的理论期望
    sim_ST_mean = np.mean(ST)
    cov = np.cov(payoff, ST)[0, 1]
    var_ST = np.var(ST)
    beta = cov / var_ST
    
    # 修正后的估计
    C_CV = C_MC - beta * (sim_ST_mean - E_ST)
    # 在列表里保存本次模拟结果
    C_CV_list.append(C_CV)

# 转换为 numpy 数组
C_MC_array = np.array(C_MC_list)
C_CV_array = np.array(C_CV_list)

# ================= 绘制直方图 =================
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 标准蒙特卡洛
ax1.hist(C_MC_array, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
ax1.axvline(C_BS, color='red', linestyle='--', linewidth=2, label=f'BS 理论价格 = {C_BS:.4f}')
ax1.axvline(np.mean(C_MC_array), color='green', linestyle='-', linewidth=2, label=f'MC 均值 = {np.mean(C_MC_array):.4f}')
ax1.set_title('标准蒙特卡洛估计分布\n(1000条路径, 重复500次)')
ax1.set_xlim(5, 15)
ax1.set_xlabel('期权价格估计')
ax1.set_ylabel('频次')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 控制变量法
ax2.hist(C_CV_array, bins=30, alpha=0.7, color='lightcoral', edgecolor='black')
ax2.axvline(C_BS, color='red', linestyle='--', linewidth=2, label=f'BS 理论价格 = {C_BS:.4f}')
ax2.axvline(np.mean(C_CV_array), color='green', linestyle='-', linewidth=2, label=f'CV 均值 = {np.mean(C_CV_array):.4f}')
ax2.set_title('控制变量法估计分布\n(1000条路径, 重复500次)')
ax2.set_xlim(5, 15)
ax2.set_xlabel('期权价格估计')
ax2.set_ylabel('频次')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 输出统计信息
print(f"Black-Scholes 理论价格: {C_BS:.6f}")
print(f"标准 MC 估计均值: {np.mean(C_MC_array):.6f}, 标准差: {np.std(C_MC_array):.6f}")
print(f"控制变量法估计均值: {np.mean(C_CV_array):.6f}, 标准差: {np.std(C_CV_array):.6f}")
print(f"方差缩减比例: {(1 - (np.std(C_CV_array)/np.std(C_MC_array))**2):.2%}")

